Trajectory Analysis: 15 Days

import copy
import sys
import os.path
import os
import numpy as np

# path_to_binaryrvg   = '/data/data_geodyn/inputs/icesat2/pre_processing/traj_files_rvg'
path_to_binaryrvg   = '/data/data_geodyn/inputs/icesat2/setups'


# iisset.2018.287.bz2
years = [2018]     #,2019]
days = np.arange(0, 367)

files_list = []
count_noletter = 0
countA = 0
countB = 0

for iyear,year in enumerate(years):
    for iday,day in enumerate(days):
#         filename = 'orbit.1807001.%d.%03d' % (year, day)
        filename = 'iisset.%d.%03d' % (year, day)

#         filename = '%d.%03d' % (year, day)
        __filename = path_to_binaryrvg+ '/'+filename
        arc = '%d.%03d' % (year, day)
        if os.path.exists(__filename):
            count_noletter+=1
            files_list.append(arc)
        if os.path.exists(__filename + 'A'):
            countA += 1
            pass
        if os.path.exists(__filename + 'B'):
            countB += 1
            pass
        if os.path.exists(__filename + '.gz'):
            count_noletter+=1
            files_list.append(arc + '')
        if os.path.exists(__filename + 'A.gz'):
            countA += 1
            pass
        if os.path.exists(__filename + 'B.gz'):
            countB += 1
            pass
        if os.path.exists(__filename + '.bz2'):
            count_noletter+=1
            files_list.append(arc + '')
        if os.path.exists(__filename + 'A.bz2'):
            countA += 1
            pass
        if os.path.exists(__filename + 'B.bz2'):
            countB += 1
            pass


arc_list = files_list[1:]

arc_list = arc_list[1:39]

print('number of arcs:', len(arc_list))
print(arc_list)

number of arcs: 38
['2018.291', '2018.292', '2018.293', '2018.294', '2018.295', '2018.296', '2018.297', '2018.298', '2018.299', '2018.300', '2018.303', '2018.304', '2018.305', '2018.306', '2018.307', '2018.308', '2018.309', '2018.312', '2018.313', '2018.314', '2018.315', '2018.316', '2018.317', '2018.318', '2018.319', '2018.320', '2018.321', '2018.322', '2018.323', '2018.324', '2018.325', '2018.326', '2018.327', '2018.328', '2018.330', '2018.331', '2018.334', '2018.335']

arc_list = arc_list[18:33]
arc_list

['2018.313',
 '2018.314',
 '2018.315',
 '2018.316',
 '2018.317',
 '2018.318',
 '2018.319',
 '2018.320',
 '2018.321',
 '2018.322',
 '2018.323',
 '2018.324',
 '2018.325',
 '2018.326',
 '2018.327']

Run GEODYN over desired arcs for all dens:

### Identify which arcs you want to run:

#------ A dictionary containing the run parameters ------
run_params = {}
run_params['arc']              =  arc_list
run_params['satellite']        =  'icesat2'
run_params['SpecialRun_name']  =  '_TrajAnalysis' #'_TrajAnalysis_fixed_cd'
run_params['verbose']          =  False
run_params['accels']           =  False
run_params['request_data']      = ['AdjustedParams',
                                    'Trajectory_orbfil',
                                    'Density',
                                    'Residuals_obs',
                                    'Residuals_summary',
                                  ]

# %load_ext autoreload
# %autoreload 2

# import pickle
# sys.path.insert(0, '/data/geodyn_proj/pygeodyn/pygeodyn_develop/')
# from PYGEODYN import Pygeodyn

# # Obj_Geodyn = {}

# for imodel,val_model in enumerate( ['msis2']):#,'msis00' ,'msis86','dtm87','jaachi71']):  # 'msis2' ,'msis86', 'msis00','dtm87',
#     run_params1 = copy.deepcopy(run_params)
#     run_params1['den_model']  =  val_model
#     run_params1['action']           =  'run'

#     ### Load the data into an object
#     Obj_Geodyn = Pygeodyn(run_params1)
#     Obj_Geodyn.RUN_GEODYN()

# import sys
# sys.exit()

Get Data

Write DataObject to a pickle

# %load_ext autoreload
# %autoreload 2

# import pickle
# sys.path.insert(0, '/data/geodyn_proj/pygeodyn/pygeodyn_develop/')
# from PYGEODYN import Pygeodyn

# # Obj_Geodyn = {}

# for imodel,val_model in enumerate( ['msis2']):#,'msis00' ,'msis86']):  # 'msis2' ,'msis86', 'msis00','dtm87', 'jaachia71'
#     read_params = copy.deepcopy(run_params)
#     read_params['den_model']  =  val_model
#     read_params['action']  =  'read'

#     ### Load the data into an object
#     Obj_Geodyn = Pygeodyn(read_params)
#     Obj_Geodyn.getData()


#     #### Pickle the object to save it
#     print('Saving pickle')
#     dir_save = '/data/data_geodyn/results/icesat2/'+'/'
#     filehandler = open(dir_save+'icesat2'+run_params['SpecialRun_name']+'_313_327_'+read_params['den_model']+'.pkl', 'wb')
#     pickle.dump(Obj_Geodyn, filehandler)
#     filehandler.close()
#     Obj_Geodyn = 0
#     print('Saved pickle')



# import sys
# sys.exit(0)

Load Pickles

# %reset out
# read the object in from a pickle
import pickle
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/pygeodyn_develop/')
from PYGEODYN import Pygeodyn

sys.path.insert(0,'/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/util_dir/')
from common_functions          import Pygeodyn_OBJECT_freeupmemory



dir_save = '/data/data_geodyn/results/icesat2/'


# Obj_Geodyn = {}
# for imodel,val_model in enumerate( ['msis86', 'msis00','dtm87','jaachia71']):
# #     read_params = copy.deepcopy(run_params)
# #     read_params['den_model']  =  val_model


#     filehandler = open(dir_save+'icesat2'+run_params['SpecialRun_name']+'_313_327_'+val_model+'.pkl', 'rb')
#     Obj_Geodyn[val_model] = pickle.load(filehandler)

#     filehandler.close()


#     print(val_model)
# # Obj_data_pkl['msis2'].__dict__.keys()

# Obj_Geodyn = Obj_data_pkl


#
############    'msis86', 'msis00','dtm87','jaachia71'

filehandler = open(dir_save+'icesat2'+'_TrajAnalysis'+'_313_327_'+'msis86'+'.pkl', 'rb')
Obj_Geodyn_msis86 = pickle.load(filehandler)
filehandler.close()
Obj_Geodyn_msis86 = Pygeodyn_OBJECT_freeupmemory(Obj_Geodyn_msis86)
print('msis86')

filehandler = open(dir_save+'icesat2'+'_TrajAnalysis'+'_313_327_'+'msis00'+'.pkl', 'rb')
Obj_Geodyn_msis00 = pickle.load(filehandler)
filehandler.close()
Obj_Geodyn_msis00 = Pygeodyn_OBJECT_freeupmemory(Obj_Geodyn_msis00)
print('msis00')

filehandler = open(dir_save+'icesat2'+'_TrajAnalysis'+'_313_327_'+'dtm87'+'.pkl', 'rb')
Obj_Geodyn_dtm87 = pickle.load(filehandler)
filehandler.close()
Obj_Geodyn_dtm87 = Pygeodyn_OBJECT_freeupmemory(Obj_Geodyn_dtm87)
print('dtm87')


filehandler = open(dir_save+'icesat2'+run_params['SpecialRun_name']+'_313_327_'+'msis2'+'.pkl', 'rb')
Obj_Geodyn_msis2 = pickle.load(filehandler)
filehandler.close()
Obj_Geodyn_msis2 = Pygeodyn_OBJECT_freeupmemory(Obj_Geodyn_msis2)
print('msis2')


filehandler = open(dir_save+'icesat2'+run_params['SpecialRun_name']+'_313_327_'+'jaachia71'+'.pkl', 'rb')
Obj_Geodyn_jaachia71 = pickle.load(filehandler)
filehandler.close()
Obj_Geodyn_jaachia71 = Pygeodyn_OBJECT_freeupmemory(Obj_Geodyn_jaachia71)
print('jaachia71')

msis86
msis00
dtm87
msis2
jaachia71

Read Memory usage:

#!/usr/bin/env python
import psutil
# gives a single float value
print(psutil.cpu_percent())
# gives an object with many fields
print(psutil.virtual_memory())
# you can convert that object to a dictionary
print(dict(psutil.virtual_memory()._asdict()))
# you can have the percentage of used RAM
print(psutil.virtual_memory().percent)
# 79.2
# you can calculate percentage of available memory
print(psutil.virtual_memory().available * 100 / psutil.virtual_memory().total)
# 20.8

0.0
svmem(total=3971731456, available=1025040384, percent=74.2, used=2707906560, free=315379712, active=2617016320, inactive=811991040, buffers=0, cached=948445184, shared=17334272, slab=62181376)
{'total': 3971731456, 'available': 1025040384, 'percent': 74.2, 'used': 2707906560, 'free': 315379712, 'active': 2617016320, 'inactive': 811991040, 'buffers': 0, 'cached': 948445184, 'shared': 17334272, 'slab': 62181376}
74.2
25.808401080377575

import gc
gc.collect()

15

Plots

import plotly.graph_objects as go
from plotly.offline import plot, iplot
from plotly.subplots import make_subplots
import plotly.express as px


config = dict({
                'displayModeBar': True,
                'responsive': False,
                'staticPlot': False,
                'displaylogo': False,
                'showTips': False,
                })

Plot: Residual Summaries

Resids = PCE - POD Trajectory

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_residual_meas_summary
from PYGEODYNAnalysis_icesat2PCEtrajectory import rms_summary_table


Obj_list = [Obj_Geodyn_msis86,
            Obj_Geodyn_msis00,
            Obj_Geodyn_msis2,
            Obj_Geodyn_dtm87,
            Obj_Geodyn_jaachia71]

rms_summary_table(Obj_list)


fig = make_subplots(rows=2, cols=1,
     subplot_titles=(["Mean Residuals per Arc", 'RMS of Fit per Arc']),
     vertical_spacing = 0.1)
fig = plot_residual_meas_summary(fig, Obj_Geodyn_msis86   , 0)
fig = plot_residual_meas_summary(fig, Obj_Geodyn_msis00   , 1)
fig = plot_residual_meas_summary(fig, Obj_Geodyn_msis2    , 2)
fig = plot_residual_meas_summary(fig, Obj_Geodyn_dtm87    , 3)
fig = plot_residual_meas_summary(fig, Obj_Geodyn_jaachia71, 4)

fig.show(config=config)

+———————————————————+—————————————————————————+————————————————+
|                     Summary Across all Arcs                  |
+———————————————————+—————————————————————————+————————————————+
+   Density Model   +   Mean Residual (cm)    +   RMS of Fit   +
+-------------------+-------------------------+----------------+
+       msis86      +       2.25333e+01       +  1.87513e-01   +
+       msis00      +       1.09333e+01       +  1.75633e-01   +
+       msis2       +       -1.04000e+01      +  1.37793e-01   +
+       dtm87       +       -2.29333e+01      +  3.89207e-01   +
+     jaachia71     +       -1.21333e+01      +  2.51780e-01   +
+———————————————————+—————————————————————————+————————————————+

Plot: XYZ Component Residuals

%load_ext autoreload
%autoreload 2

from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_residuals_observed

fig = make_subplots(rows=3, cols=1,
            subplot_titles=(['X', 'Y', 'Z']),
            vertical_spacing = 0.1,
                       )

fig = plot_residuals_observed(fig, Obj_Geodyn_msis86    , 0)
fig = plot_residuals_observed(fig, Obj_Geodyn_msis00    , 1)
fig = plot_residuals_observed(fig, Obj_Geodyn_msis2     , 2)
fig = plot_residuals_observed(fig, Obj_Geodyn_dtm87     , 3)
fig = plot_residuals_observed(fig, Obj_Geodyn_jaachia71 , 4)

fig.update_layout(title="Observation Residuals (PCE - Observed , T.O.R.)")

# fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Plot: Adjusted Cd

%load_ext autoreload
%autoreload 2
from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_cd_and_percdiff_from_apriori


fig = make_subplots(
    rows=2, cols=1,
    subplot_titles=(["Timeseries of Cd", "Ratio of Adjusted Cd to a priori (Cd=2.2)"]),
    vertical_spacing = 0.08,
    )
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn_msis86     , 0)
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn_msis00     , 1)
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn_msis2      , 2)
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn_dtm87      , 3)
fig = plot_cd_and_percdiff_from_apriori(fig,  Obj_Geodyn_jaachia71  , 4)


fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

obj_m1  = Obj_Geodyn_msis86
plot_num = 0
import sys
sys.path.insert(0,'/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/util_dir/')
from common_functions          import Convert_cartesian_to_RSW, Convert_cartesian_to_NTW
import pandas as pd


###### GET THE PCE DATA:
# StateVector_PCE_datafile = '/data/data_geodyn/inputs/icesat2/setups/PCE_ascii.txt'
#     os.system('bunzip2 -v '+ StateVector_epochs_datafile +'.bz2')
import plotly.graph_objects as go
from plotly.offline import plot, iplot
from plotly.subplots import make_subplots
import plotly.express as px


config = dict({
                'displayModeBar': True,
                'responsive': False,
                'staticPlot': True,
                'displaylogo': False,
                'showTips': False,
                })





from datetime import datetime,timedelta

# Simplify Plotting Schemes:
col1 = px.colors.qualitative.Plotly[0]
col2 = px.colors.qualitative.Plotly[1]
col3 = px.colors.qualitative.Plotly[2]
col4 = px.colors.qualitative.Plotly[3]
col5 = px.colors.qualitative.Plotly[4]
col6 = px.colors.qualitative.Plotly[5]

# Obj_Geodyn_msis86.__dict__['Trajectory_orbfil']['2018.327']['data_record']

# last_nonpredicted_time = pd.Series(pd.to_datetime(obj_m1.__dict__['Residuals_obs'][first_arc]['Date'].iloc[-1]))

# dates = obj_m1.__dict__['Trajectory_orbfil'][first_arc]['data_record']['Date_UTC']


# (dates - last_nonpredicted_time).abs().min()

# (dates - last_nonpredicted_time)
# print(dates.index.get_loc(last_nonpredicted_time, method='nearest'))
# dates.index[dates.index.get_loc(last_nonpredicted_time, method='nearest')]



# dates.index.get_loc(last_nonpredicted_time, method='nearest')
# dates.iloc[dates.index.get_loc(last_nonpredicted_time[0], method='nearest')]
# last_nonpredicted_time
# last_index = int(np.round(np.size(dates)/))
# dates[:17821]

# 195000*15

Plot: XYZ+NTW+RSW

# %load_ext autoreload
# %autoreload 2
# # def ARCOVERLAP_2arcs_ObsResids_NTW_intrack(fig, obj_m1, plot_num, arc1, arc2):

# from common_functions          import Convert_cartesian_to_NTW_returnall, Convert_cartesian_to_RSW_returnall


# fig = make_subplots(rows=3, cols=3,
#             shared_xaxes=True,
# #             subplot_titles=(('X','N','R', 'Y','T','S', 'Z','W', 'W')),
# #             subplot_titles=([ ['X', 'Y','Z'],
# #                             ['N', 'T','W']]),
#             column_titles = (['XYZ', 'NTW','RSW']),
#             vertical_spacing = 0.081,
#             horizontal_spacing = 0.081,

#                    )

# #### ------------------------ plotting stuff
# model_m1 = obj_m1.__dict__['global_params']['den_model']
# if plot_num == 0:
#     col = col1
#     x_annot = 1.05
#     y_annot = .97
#     m_size = 3
# elif plot_num == 1:
#     x_annot = 1.05
#     y_annot = .8
#     col = col2
#     m_size = 2.5
# elif plot_num == 2:
#     x_annot = 1.05
#     y_annot = .55
#     col = col3
#     m_size = 2
# data_skip = 14
# #### ------------------------ plotting stuff


# #### ----------------------------------------------------------------------
# #####  read the PCE data but only for the lines during which we have data:
# #### the below block of code just READS the PCE ascii at the times in the ORBFIL
# # first_arc = obj_m1.__dict__['global_params']['arc_input'][0]
# # last_arc  = obj_m1.__dict__['global_params']['arc_input'][-1]


# # first_arc_first_time = obj_m1.__dict__['Trajectory_orbfil'][first_arc]['data_record']['Date_UTC'].iloc[0],
# # last_arc_last_time   = obj_m1.__dict__['Trajectory_orbfil'][first_arc]['data_record']['Date_UTC'].iloc[-1],
# # first_arc_first_time_str =  str(first_arc_first_time[0])
# # last_arc_last_time_str   =  str(  last_arc_last_time[0])

# # with open(StateVector_PCE_datafile, 'r') as f:
# #     for line_no, line_text in enumerate(f):
# #         if first_arc_first_time_str in line_text:
# #             print('first ',line_text, line_no)
# #             first_line = line_no
# #         if last_arc_last_time_str in line_text:
# #             print('last ',line_text, line_no)
# #             last_line = line_no
# #             break
# # ### load PCE data as a dataframe
# # PCE_data = pd.read_csv(StateVector_PCE_datafile,
# #             skiprows = first_line,
# #             nrows=last_line- first_line,
# # #            nrows=2925000,
# #             sep = '\s+',
# #             dtype=str,
# #             names = [
# #                     'Date',
# #                     'MJDSECs',
# #                     'RSECS', #(fractional secs)
# #                     'GPS offset', # (UTC - GPS offset (secs))
# #                     'X',
# #                     'Y',
# #                     'Z',
# #                     'X_dot',
# #                     'Y_dot',
# #                     'Z_dot',
# #                     'YYMMDDhhmmss',
# #                         ],)

# # print('data loaded')

# # #     os.system('bzip2 -v '+ StateVector_epochs_datafile )
# # #### ----------------------------------------------------------------------



# # ####---------------------------------------------------------
# # ### Collect the PCE data (true) and the Orbfil data (estimated)
# # #### use pd.merge to only assess them where they have common dates.
# # PCE_data['Date_pd'] = pd.to_datetime(PCE_data['Date'])


# for iarc,arcval in enumerate([obj_m1.__dict__['global_params']['arc_input'][0]]):
#     arc_first_time = obj_m1.__dict__['Trajectory_orbfil'][arcval]['data_record']['Date_UTC'].iloc[0],
#     arc_last_time   = obj_m1.__dict__['Trajectory_orbfil'][arcval]['data_record']['Date_UTC'].iloc[-1],
#     arc_first_time_str =  str(arc_first_time[0])
#     arc_last_time_str   =  str(  arc_last_time[0])

#     with open(StateVector_PCE_datafile, 'r') as f:
#         for line_no, line_text in enumerate(f):
#             if arc_first_time_str in line_text:
#                 print('first ',line_text, line_no)
#                 first_line = line_no
#             if arc_last_time_str in line_text:
#                 print('last ',line_text, line_no)
#                 last_line = line_no
#                 break
#     ### load PCE data as a dataframe
#     PCE_data = pd.read_csv(StateVector_PCE_datafile,
#             skiprows = first_line,
#             nrows=last_line- first_line,
#     #            nrows=2925000,
#             sep = '\s+',
#             dtype=str,
#             names = [
#                     'Date',
#                     'MJDSECs',
#                     'RSECS', #(fractional secs)
#                     'GPS offset', # (UTC - GPS offset (secs))
#                     'X',
#                     'Y',
#                     'Z',
#                     'X_dot',
#                     'Y_dot',
#                     'Z_dot',
#                     'YYMMDDhhmmss',
#                         ],)

#     print('data loaded')

#     #     os.system('bzip2 -v '+ StateVector_epochs_datafile )
#     #### ----------------------------------------------------------------------



#     ####---------------------------------------------------------
#     ### Collect the PCE data (true) and the Orbfil data (estimated)
#     #### use pd.merge to only assess them where they have common dates.
#     PCE_data['Date_pd'] = pd.to_datetime(PCE_data['Date'])



#     ### load orbfil for this arc as a dataframe
#     orbfil_arc = obj_m1.__dict__['Trajectory_orbfil'][arcval]['data_record']
#     orbfil_arc['Date_pd'] = pd.to_datetime(orbfil_arc ['Date_UTC'])


#     C_1 = pd.merge(left=orbfil_arc, left_on='Date_pd',
#          right=PCE_data, right_on='Date_pd')

#     print('data merged')

#     del PCE_data


#     print('data deleted')
#     X    = C_1['Satellite Inertial X coordinate']
#     Y    = C_1['Satellite Inertial Y coordinate']
#     Z    = C_1['Satellite Inertial Z coordinate']
#     Xdot = C_1['Satellite Inertial X velocity']
#     Ydot = C_1['Satellite Inertial Y velocity']
#     Zdot = C_1['Satellite Inertial Z velocity']

#     #####======================================================================
#     ##### PLOT THE XYZ COMPONENTS!!!
#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][:17821:data_skip],
#                                y=X[:17821:data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=1, col=1,
#                              )
#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][:17821:data_skip],
#                                y=Y[:17821:data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=2, col=1,
#                              )
#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][:17821:data_skip],
#                                y=Z[:17821:data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=3, col=1,
#                              )
#     #####======================================================================
#     ##### PLOT THE NTW COMPONENTS!!!
#     state_vector = np.transpose(np.array([X, Y, Z, Xdot, Ydot, Zdot]))
#     ntw_comps = [Convert_cartesian_to_NTW_returnall(x) for x in state_vector]

#     N = np.transpose(ntw_comps)[0]
#     T = np.transpose(ntw_comps)[1]
#     W = np.transpose(ntw_comps)[2]

#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][:17821:data_skip],
#                                y=N[:17821:data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=1, col=2,
#                              )
#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][:17821:data_skip],
#                                y=T[:17821:data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=2, col=2,
#                              )
#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][:17821:data_skip],
#                                y=W[:17821:data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=3, col=2,
#                              )

#     #####======================================================================
#     ##### PLOT THE RSW COMPONENTS!!!
#     state_vector = np.transpose(np.array([X, Y, Z, Xdot, Ydot, Zdot]))
#     rsw_comps = [Convert_cartesian_to_RSW_returnall(x) for x in state_vector]

#     R = np.transpose(rsw_comps)[0]
#     S = np.transpose(rsw_comps)[1]
#     W = np.transpose(rsw_comps)[2]

#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][:17821:data_skip],
#                                y=R[:17821:data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=1, col=3,
#                              )
#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][:17821:data_skip],
#                                y=S[:17821:data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=2, col=3,
#                              )
#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][:17821:data_skip],
#                                y=W[:17821:data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=3, col=3,
#                              )



# # ### Start of second arc
# # overlap_start = obj_m1.__dict__['Residuals_obs'][arc1]['Date'].iloc[-1]
# # ### End of first arc
# # overlap_end   = obj_m1.__dict__['Trajectory_orbfil'][arc1]['data_record']['Date_UTC'].iloc[-1]
# # fig.add_vrect(  x0=overlap_start, x1=overlap_end,
# #                 fillcolor='LightSkyBlue', opacity=0.2,
# #                 layer="below", line_width=0)


# fig.update_layout(title="Coordinate System Comparisons (all in meters)")
# fig.update_layout(
#                 autosize=False,
#                 width=1500,
#                 height=800,
#                 font=dict(size=12),
#                 legend= {'itemsizing': 'constant'})

# fig.update_yaxes( title="X",exponentformat= 'power',row=1, col=1)
# fig.update_yaxes( title="Y",exponentformat= 'power',row=2, col=1)
# fig.update_yaxes( title="Z",exponentformat= 'power',row=3, col=1)

# fig.update_yaxes( title="N",exponentformat= 'power',row=1, col=2)
# fig.update_yaxes( title="T",exponentformat= 'power',row=2, col=2)
# fig.update_yaxes( title="W",exponentformat= 'power',row=3, col=2)

# fig.update_yaxes( title="R",exponentformat= 'power',row=1, col=3)
# fig.update_yaxes( title="S",exponentformat= 'power',row=2, col=3)
# fig.update_yaxes( title="W",exponentformat= 'power',row=3, col=3)



# fig.update_xaxes( title="Date", row=3, col=1)
# fig.update_xaxes( title="Date", row=3, col=2)
# fig.update_xaxes( title="Date", row=3, col=3)
# # fig.update_yaxes(title_text="Residuals (cm)", row=1, col=1, secondary_y=True, color='SkyBlue')

Plot: NTW Residuals w/ CD Scaling factor

# #!/usr/bin/env python
# import psutil
# # gives a single float value
# print(psutil.cpu_percent())
# # gives an object with many fields
# print(psutil.virtual_memory())
# # you can convert that object to a dictionary
# print(dict(psutil.virtual_memory()._asdict()))
# # you can have the percentage of used RAM
# print(psutil.virtual_memory().percent)
# # 79.2
# # you can calculate percentage of available memory
# print(psutil.virtual_memory().available * 100 / psutil.virtual_memory().total)
# # 20.8

%load_ext autoreload
%autoreload 2
from PYGEODYNAnalysis_icesat2PCEtrajectory import NTW_CDratio_IntrackResids


fig = make_subplots(rows=2, cols=1,
            shared_xaxes=True,
            subplot_titles=(['Adjusted Cd ratio to a priori (2.2)', 'In-Track Component Residuals (PCE-ORBFIL)']),
            vertical_spacing = 0.1,
            specs=[ [{"secondary_y": False }],
                    [{"secondary_y": False }]],)


fig = NTW_CDratio_IntrackResids(fig,  Obj_Geodyn_msis86     , 0)
fig = NTW_CDratio_IntrackResids(fig,  Obj_Geodyn_msis00     , 1)
fig = NTW_CDratio_IntrackResids(fig,  Obj_Geodyn_msis2      , 2)
fig = NTW_CDratio_IntrackResids(fig,  Obj_Geodyn_dtm87      , 3)
fig = NTW_CDratio_IntrackResids(fig,  Obj_Geodyn_jaachia71  , 4)


fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
2018.313
2018.315
2018.317
2018.319
2018.321
2018.323
2018.325
2018.327
2018.313
2018.315
2018.317
2018.319
2018.321
2018.323
2018.325
2018.327
2018.313
2018.315
2018.317
2018.319
2018.321
2018.323
2018.325
2018.327
2018.313
2018.315
2018.317
2018.319
2018.321
2018.323
2018.325
2018.327
2018.313
2018.315
2018.317
2018.319
2018.321
2018.323
2018.325
2018.327

In-track arc overlap

def NTW_CDratio_IntrackResids_arcoverlap(fig, obj_m1, plot_num):

#     import sys
#     sys.path.insert(0,'/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/util_dir/')
#     from common_functions          import Convert_cartesian_to_RSW, Convert_cartesian_to_NTW


#     fig = make_subplots(rows=2, cols=1,
#                 shared_xaxes=True,
#                 subplot_titles=(['Adjusted Cd ratio to a priori (2.2)', 'In-Track Component Residuals (PCE-ORBFIL)']),
#                 vertical_spacing = 0.1,
#                 specs=[ [{"secondary_y": False }],
#                         [{"secondary_y": False }]],)

    model_m1 = obj_m1.__dict__['global_params']['den_model']

    if plot_num == 0:
        col = col1
        x_annot = 1.09
        y_annot = .90
        m_size = 5
    elif plot_num == 1:
        x_annot = 1.09
        y_annot = .7
        col = col2
        m_size = 5
    elif plot_num == 2:
        x_annot = 1.09
        y_annot = .5
        col = col3
        m_size = 5
    elif plot_num == 3:
        x_annot = 1.09
        y_annot = .3
        col = col4
        m_size = 5
    elif plot_num == 4:
        x_annot = 1.09
        y_annot = .1
        col = col5
        m_size = 5

    ###### GET THE PCE DATA:
    StateVector_PCE_datafile = '/data/data_geodyn/inputs/icesat2/setups/PCE_ascii.txt'
    SAT_ID = int(obj_m1.__dict__['global_params']['SATID'])
    which_stat = 'CURRENT_VALUE'
    data_skip = 15

    ####--------------------- Residual  ---------------------
    for ii,arc in enumerate(obj_m1.__dict__['global_params']['arc_input'][:2]):
        print(arc)

        if ii == 1:
            mark_symbol = 'triangle-up'
        else:
            mark_symbol = 'circle'

    #     arc_first_time = obj_m1.__dict__['Trajectory_orbfil'][arcval]['data_record']['Date_UTC'].iloc[0],
    #     arc_last_time   = obj_m1.__dict__['Trajectory_orbfil'][arcval]['data_record']['Date_UTC'].iloc[-1],
    #     arc_first_time_str =  str(arc_first_time[0])
    #     arc_last_time_str   =  str(  arc_last_time[0])

        arc_first_time  = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']['Date_UTC'].iloc[0],
        arc_last_time   = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']['Date_UTC'].iloc[-1]
        arc_first_time_str     =  str(arc_first_time[0])#.replace( "'",' ')
        arc_last_time_str      =  str(arc_last_time)#.replace( "'",' ')

        # print('first_arc_first_time',first_arc_first_time)
        ####---------------------------------------------------------
        with open(StateVector_PCE_datafile, 'r') as f:
            for line_no, line_text in enumerate(f):

                if arc_first_time_str in line_text:
                    first_line = line_no
    #                 print('first worked', arc_first_time_str)
                elif arc_last_time_str in line_text:
                    last_line = line_no
    #                 print('last worked', arc_last_time_str)

                    break
        PCE_data = pd.read_csv(StateVector_PCE_datafile,
                    skiprows = first_line,
                    nrows=last_line-first_line,
                    sep = '\s+',
                    dtype=str,
                    names = [
                            'Date',
                            'MJDSECs',
                            'RSECS', #(fractional secs)
                            'GPS offset', # (UTC - GPS offset (secs))
                            'X',
                            'Y',
                            'Z',
                            'X_dot',
                            'Y_dot',
                            'Z_dot',
                            'YYMMDDhhmmss',
                                ],)

        PCE_data['Date_pd'] = pd.to_datetime(PCE_data['Date'])
        orbfil_arc1 = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']
        orbfil_arc1['Date_pd'] = pd.to_datetime(orbfil_arc1 ['Date_UTC'])

        C_1 = pd.merge(left=orbfil_arc1, left_on='Date_pd',
             right=PCE_data, right_on='Date_pd')

        X = C_1['Satellite Inertial X coordinate']
        Y = C_1['Satellite Inertial Y coordinate']
        Z = C_1['Satellite Inertial Z coordinate']
        Xdot = C_1['Satellite Inertial X velocity']
        Ydot = C_1['Satellite Inertial Y velocity']
        Zdot = C_1['Satellite Inertial Z velocity']
        state_vector = np.transpose(np.array([X, Y, Z, Xdot, Ydot, Zdot]))
        InTrack_comp_orbfil = [Convert_cartesian_to_NTW(x) for x in state_vector]

        X = C_1['X'].astype(float)
        Y = C_1['Y'].astype(float)
        Z = C_1['Z'].astype(float)
        Xdot = C_1['X_dot'].astype(float)
        Ydot = C_1['Y_dot'].astype(float)
        Zdot = C_1['Z_dot'].astype(float)
        state_vector = np.transpose(np.array([X, Y, Z, Xdot, Ydot, Zdot]))
        InTrack_comp_PCE = [Convert_cartesian_to_NTW(x) for x in state_vector]

        resid = (np.array(InTrack_comp_PCE) - np.array(InTrack_comp_orbfil))*1e2

    #     print('C_1              ', np.size(C_1['Date_pd']))
    #     print('resid            ', np.size(resid))
    #     print('InTrack_comp_PCE ', np.size(InTrack_comp_PCE))

        fig.add_trace(go.Scattergl(x=C_1['Date_pd'][::data_skip],
                                   y=resid[::data_skip],
                                 name= '(PCE-orbfil)',
                                 mode='markers',
                                 marker=dict(color=col,
                                 size=m_size, symbol= mark_symbol,),
                                 showlegend=False,
                                 ),
                                 secondary_y=False,
                                 row=2, col=1,
                                 )

        str_run_param = 'run_parameters'+ arc
        final_iter = obj_m1.__dict__[str_run_param]['str_iteration']

        i_arc = ii+1
        last_iter = list(obj_m1.AdjustedParams[arc].keys())[-1]
        time_dep_cd_dates = list(obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'].keys())


        val_list_1 = []
        for i in time_dep_cd_dates:
    #         print(obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
            val_list_1.append(obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
    #         print('CD_DATE  ',pd.to_datetime(i))
            cd_ratio =  obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'][i][which_stat]/ 2.2
            fig.add_trace(go.Scattergl(x=  [pd.to_datetime(i)],
                                       y=  [cd_ratio],
                               name= model_m1,
#                                mode='markers+text',
                               mode='markers',
                               marker=dict(
                               symbol= mark_symbol,
                               color=col,
                               size=7,
                               ),
#                                text =  ''+str(round(obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'][i][which_stat], 2)), #str(cd_ratio),
#                                textposition="top center",
                               showlegend=False,
                               ),
                               row=1, col=1,)

        add_dt = pd.to_timedelta(180,'m')
        overlap_end   = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']['Date_UTC'].iloc[-1]

        date = pd.to_datetime(i)

#         while date < overlap_end:
#             date += add_dt
#             fig.add_trace(go.Scattergl(x= [date] ,
#                                        y=  list(np.ones(1)),
#                                        name= model_m1,
#                                        mode='markers+text',
#                                        marker=dict(
#                                        color=col,
#                                        size=7,
#                                        ),
#                                        text =  '2.2', #str(cd_ratio),
#                                        textposition="bottom center",
#                                        showlegend=False,
#                                        ),
#                                        row=1, col=1,)


        ### Start of second arc
        overlap_start = obj_m1.__dict__['Residuals_obs'][arc]['Date'].iloc[-1]
        ### End of first arc
        overlap_end   = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']['Date_UTC'].iloc[-1]
        fig.add_vrect(  x0=overlap_start, x1=overlap_end,
                        fillcolor='LightSkyBlue', opacity=0.2,
                        layer="below", line_width=0)

    # fig = legend_as_annotation(fig, obj_m1.__dict__['global_params']['den_model'], col, x_annot, y_annot)

    fig.update_layout(title="NTW Coord. System + Predicted Window (light blue window)")
    fig.update_layout(
                    autosize=False,
                    width=900,
                    height=700,
                    font=dict(size=12),
                    legend= {'itemsizing': 'constant'})

    fig.update_yaxes( title="Cd Ratio (Adjusted Cd/2.2) ",exponentformat= 'power',row=1, col=1)
    fig.update_yaxes( title="Residual (cm)",       exponentformat= 'power',row=2, col=1)
    fig.update_xaxes( title="Date", row=2, col=1)
    fig.update_yaxes(title_text="Residuals (cm)", row=1, col=1, secondary_y=True, color='SkyBlue')
    fig.update_traces(textfont_size=10, textfont_color=col)
#     fig.update_xaxes(range=[obj_m1.__dict__['Trajectory_orbfil'][obj_m1.__dict__['global_params']['arc_input'][0]]['data_record']['Date_UTC'].iloc[0], C_1['Date_pd'].iloc[-1]+ pd.to_timedelta(60,'m')])

    return(fig)


%load_ext autoreload
%autoreload 2
# from PYGEODYNAnalysis_icesat2PCEtrajectory import NTW_CDratio_IntrackResids_arcoverlap


fig = make_subplots(rows=2, cols=1,
            shared_xaxes=True,
            subplot_titles=(['Adjusted Cd ratio to a priori (2.2)', 'In-Track Component Residuals (PCE-ORBFIL)']),
            vertical_spacing = 0.1,
            specs=[ [{"secondary_y": False }],
                    [{"secondary_y": False }]],)


fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_msis86     , 0)
# fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_msis00     , 1)
# fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_msis2      , 2)
fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_dtm87      , 3)
fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_jaachia71  , 4)


fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
2018.313
2018.314
2018.313
2018.314
2018.313
2018.314

%load_ext autoreload
%autoreload 2
# from PYGEODYNAnalysis_icesat2PCEtrajectory import NTW_CDratio_IntrackResids_arcoverlap


fig = make_subplots(rows=2, cols=1,
            shared_xaxes=True,
            subplot_titles=(['Adjusted Cd ratio to a priori (2.2)', 'In-Track Component Residuals (PCE-ORBFIL)']),
            vertical_spacing = 0.1,
            specs=[ [{"secondary_y": False }],
                    [{"secondary_y": False }]],)


fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_msis86     , 0)
fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_msis00     , 1)
fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_msis2      , 2)
# fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_dtm87      , 3)
# fig = NTW_CDratio_IntrackResids_arcoverlap(fig,  Obj_Geodyn_jaachia71  , 4)


fig.show(config=config)

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
2018.313
2018.314
2018.313
2018.314
2018.313
2018.314

# # def ARCOVERLAP_2arcs_ObsResids_NTW_intrack(fig, obj_m1, plot_num, arc1, arc2):

# import sys
# sys.path.insert(0,'/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/util_dir/')
# from common_functions          import Convert_cartesian_to_RSW, Convert_cartesian_to_NTW


# fig = make_subplots(rows=2, cols=1,
#             shared_xaxes=True,
#             subplot_titles=(['Adjusted Cd ratio to a priori (2.2)', 'In-Track Component Residuals (PCE-ORBFIL)']),
#             vertical_spacing = 0.1,
#             specs=[ [{"secondary_y": False }],
#                     [{"secondary_y": False }]],)

# model_m1 = obj_m1.__dict__['global_params']['den_model']

# if plot_num == 0:
#     col = col1
#     x_annot = 1.05
#     y_annot = .97
#     m_size = 3
# elif plot_num == 1:
#     x_annot = 1.05
#     y_annot = .8
#     col = col2
#     m_size = 2.5
# elif plot_num == 2:
#     x_annot = 1.05
#     y_annot = .55
#     col = col3
#     m_size = 2

# ###### GET THE PCE DATA:
# StateVector_PCE_datafile = '/data/data_geodyn/inputs/icesat2/setups/PCE_ascii.txt'
# SAT_ID = int(obj_m1.__dict__['global_params']['SATID'])
# which_stat = 'CURRENT_VALUE'
# data_skip = 15

# ####--------------------- Residual  ---------------------
# for ii,arc in enumerate(obj_m1.__dict__['global_params']['arc_input'][::2]):
#     print(arc)



# #     arc_first_time = obj_m1.__dict__['Trajectory_orbfil'][arcval]['data_record']['Date_UTC'].iloc[0],
# #     arc_last_time   = obj_m1.__dict__['Trajectory_orbfil'][arcval]['data_record']['Date_UTC'].iloc[-1],
# #     arc_first_time_str =  str(arc_first_time[0])
# #     arc_last_time_str   =  str(  arc_last_time[0])

#     arc_first_time  = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']['Date_UTC'].iloc[0],
#     arc_last_time   = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']['Date_UTC'].iloc[-1]
#     arc_first_time_str     =  str(arc_first_time[0])#.replace( "'",' ')
#     arc_last_time_str      =  str(arc_last_time)#.replace( "'",' ')

#     # print('first_arc_first_time',first_arc_first_time)
#     ####---------------------------------------------------------
#     with open(StateVector_PCE_datafile, 'r') as f:
#         for line_no, line_text in enumerate(f):

#             if arc_first_time_str in line_text:
#                 first_line = line_no
# #                 print('first worked', arc_first_time_str)
#             elif arc_last_time_str in line_text:
#                 last_line = line_no
# #                 print('last worked', arc_last_time_str)

#                 break
#     PCE_data = pd.read_csv(StateVector_PCE_datafile,
#                 skiprows = first_line,
#                 nrows=last_line-first_line,
#                 sep = '\s+',
#                 dtype=str,
#                 names = [
#                         'Date',
#                         'MJDSECs',
#                         'RSECS', #(fractional secs)
#                         'GPS offset', # (UTC - GPS offset (secs))
#                         'X',
#                         'Y',
#                         'Z',
#                         'X_dot',
#                         'Y_dot',
#                         'Z_dot',
#                         'YYMMDDhhmmss',
#                             ],)

#     PCE_data['Date_pd'] = pd.to_datetime(PCE_data['Date'])
#     orbfil_arc1 = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']
#     orbfil_arc1['Date_pd'] = pd.to_datetime(orbfil_arc1 ['Date_UTC'])

#     C_1 = pd.merge(left=orbfil_arc1, left_on='Date_pd',
#          right=PCE_data, right_on='Date_pd')

#     X = C_1['Satellite Inertial X coordinate']
#     Y = C_1['Satellite Inertial Y coordinate']
#     Z = C_1['Satellite Inertial Z coordinate']
#     Xdot = C_1['Satellite Inertial X velocity']
#     Ydot = C_1['Satellite Inertial Y velocity']
#     Zdot = C_1['Satellite Inertial Z velocity']
#     state_vector = np.transpose(np.array([X, Y, Z, Xdot, Ydot, Zdot]))
#     InTrack_comp_orbfil = [Convert_cartesian_to_NTW(x) for x in state_vector]

#     X = C_1['X'].astype(float)
#     Y = C_1['Y'].astype(float)
#     Z = C_1['Z'].astype(float)
#     Xdot = C_1['X_dot'].astype(float)
#     Ydot = C_1['Y_dot'].astype(float)
#     Zdot = C_1['Z_dot'].astype(float)
#     state_vector = np.transpose(np.array([X, Y, Z, Xdot, Ydot, Zdot]))
#     InTrack_comp_PCE = [Convert_cartesian_to_NTW(x) for x in state_vector]

#     resid = (np.array(InTrack_comp_PCE) - np.array(InTrack_comp_orbfil))*1e2

# #     print('C_1              ', np.size(C_1['Date_pd']))
# #     print('resid            ', np.size(resid))
# #     print('InTrack_comp_PCE ', np.size(InTrack_comp_PCE))

#     fig.add_trace(go.Scattergl(x=C_1['Date_pd'][::data_skip],
#                                y=resid[::data_skip],
#                              name= '(PCE-orbfil)',
#                              mode='markers',
#                              marker=dict(color=col,
#                              size=m_size,),
#                              showlegend=False,
#                              ),
#                              secondary_y=False,
#                              row=2, col=1,
#                              )

#     str_run_param = 'run_parameters'+ arc
#     final_iter = obj_m1.__dict__[str_run_param]['str_iteration']

#     i_arc = ii+1
#     last_iter = list(obj_m1.AdjustedParams[arc].keys())[-1]
#     time_dep_cd_dates = list(obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'].keys())


#     val_list_1 = []
#     for i in time_dep_cd_dates:
# #         print(obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
#         val_list_1.append(obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
# #         print('CD_DATE  ',pd.to_datetime(i))
#         cd_ratio =  obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'][i][which_stat]/ 2.2
#         fig.add_trace(go.Scattergl(x=  [pd.to_datetime(i)],
#                                    y=  [cd_ratio],
#                            name= model_m1,
#                            mode='markers+text',
#                            marker=dict(
#                            color=col,
#                            size=7,
#                            ),
#                            text =  ''+str(round(obj_m1.AdjustedParams[arc][last_iter][SAT_ID]['0CD'][i][which_stat], 2)), #str(cd_ratio),
#                            textposition="top center",
#                            showlegend=False,
#                            ),
#                            row=1, col=1,)

#     add_dt = pd.to_timedelta(180,'m')
#     overlap_end   = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']['Date_UTC'].iloc[-1]

#     date = pd.to_datetime(i)

#     while date < overlap_end:
#         date += add_dt
#         fig.add_trace(go.Scattergl(x= [date] ,
#                                    y=  list(np.ones(1)),
#                                    name= model_m1,
#                                    mode='markers+text',
#                                    marker=dict(
#                                    color=col,
#                                    size=7,
#                                    ),
#                                    text =  '2.2', #str(cd_ratio),
#                                    textposition="bottom center",
#                                    showlegend=False,
#                                    ),
#                                    row=1, col=1,)



#     ### Start of second arc
#     overlap_start = obj_m1.__dict__['Residuals_obs'][arc]['Date'].iloc[-1]
#     ### End of first arc
#     overlap_end   = obj_m1.__dict__['Trajectory_orbfil'][arc]['data_record']['Date_UTC'].iloc[-1]
#     fig.add_vrect(  x0=overlap_start, x1=overlap_end,
#                     fillcolor='LightSkyBlue', opacity=0.2,
#                     layer="below", line_width=0)

# # fig = legend_as_annotation(fig, obj_m1.__dict__['global_params']['den_model'], col, x_annot, y_annot)

# fig.update_layout(title="NTW Coord. System + Predicted Window (light blue window)")
# fig.update_layout(
#                 autosize=False,
#                 width=900,
#                 height=700,
#                 font=dict(size=12),
#                 legend= {'itemsizing': 'constant'})

# fig.update_yaxes( title="Cd Ratio (Adjusted Cd/2.2) ",exponentformat= 'power',row=1, col=1)
# fig.update_yaxes( title="Residual (cm)",       exponentformat= 'power',row=2, col=1)
# fig.update_xaxes( title="Date", row=2, col=1)
# fig.update_yaxes(title_text="Residuals (cm)", row=1, col=1, secondary_y=True, color='SkyBlue')
# fig.update_traces(textfont_size=10, textfont_color=col)
# fig.update_xaxes(range=[obj_m1.__dict__['Trajectory_orbfil'][obj_m1.__dict__['global_params']['arc_input'][0]]['data_record']['Date_UTC'].iloc[0], C_1['Date_pd'].iloc[-1]])

# #     return(fig)

Plot: Model Density + Scaled with Cd factor

# ### Clear up some memory space...
# # def Pygeodyn_OBJECT_freeupmemory(OBJ):
# for OBJ in [Obj_Geodyn_msis86,
#             Obj_Geodyn_msis00,
#             Obj_Geodyn_dtm87,
#             Obj_Geodyn_msis2,
#             Obj_Geodyn_jaachia71]:

#         del OBJ.__dict__['Residuals_obs']
#         del OBJ.__dict__['Residuals_summary']
#         del OBJ.__dict__['Trajectory_orbfil']

# Obj_Geodyn_dtm87.__dict__.keys()

# #!/usr/bin/env python
# import psutil
# # gives a single float value
# print(psutil.cpu_percent())
# # gives an object with many fields
# print(psutil.virtual_memory())
# # you can convert that object to a dictionary
# print(dict(psutil.virtual_memory()._asdict()))
# # you can have the percentage of used RAM
# print(psutil.virtual_memory().percent)
# # 79.2
# # you can calculate percentage of available memory
# print(psutil.virtual_memory().available * 100 / psutil.virtual_memory().total)
# # 20.8

# %load_ext autoreload
# %autoreload 2

# from PYGEODYNAnalysis_icesat2PCEtrajectory import plot_ScaleDensity_with_CdScaleFactor__2

# fig = make_subplots(rows=2, cols=1,
#                 subplot_titles=(["Model Ouptut Density", "Model Density * Cd Scaling Factor"]),
#                 shared_yaxes=True,
#                 vertical_spacing = 0.1,
#                 specs=[
#                 [{"secondary_y": False}],
#                 [{"secondary_y": False}], ])


# fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn_msis86     , 0, 200)
# fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn_msis00     , 1, 200)
# fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn_msis2      , 2, 200)
# fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn_dtm87      , 3, 200)
# fig = plot_ScaleDensity_with_CdScaleFactor__2(fig,  Obj_Geodyn_jaachia71  , 4, 200)

# # min_y = 1*1e-16
# # max_y = 9*1e-12
# # fig.update_yaxes(range=[min_y, max_y], row=1, col=1)
# # fig.update_yaxes(range=[min_y, max_y], row=2, col=1)

# fig.show(config=config)

# #!/usr/bin/env python
# import psutil
# # gives a single float value
# print(psutil.cpu_percent())
# # gives an object with many fields
# print(psutil.virtual_memory())
# # you can convert that object to a dictionary
# print(dict(psutil.virtual_memory()._asdict()))
# # you can have the percentage of used RAM
# print(psutil.virtual_memory().percent)
# # 79.2
# # you can calculate percentage of available memory
# print(psutil.virtual_memory().available * 100 / psutil.virtual_memory().total)
# # 20.8

ADDITIONAL PLOTS:

We want to show the residuals in the overlap time with the PCE data subtracted away.

# ####  ARC_OVERLAP_ObsResids_XYZ

# %load_ext autoreload
# %autoreload 2

# from PYGEODYNAnalysis_icesat2PCEtrajectory import ARCOVERLAP_2arcs_ObsResids_XYZ

# fig = make_subplots(rows=3, cols=1,
#             subplot_titles=(['X', 'Y', 'Z']),
#             vertical_spacing = 0.1,
#             specs=[ [{"secondary_y": True }],
#                     [{"secondary_y": True }],
#                     [{"secondary_y": True }], ],)

# arc1 = '2018.319'  # '2018.314'
# arc2 = arc1

# fig = ARCOVERLAP_2arcs_ObsResids_XYZ(fig, Obj_Geodyn['msis2'], 0, arc1, arc2)
# # fig = ARCOVERLAP_2arcs_ObsResids_XYZ(fig, Obj_Geodyn3, 1, arc1, arc2)
# # fig = ARCOVERLAP_2arcs_ObsResids_XYZ(fig, Obj_Geodyn1, 2, arc1, arc2)

# fig.show(config=config)

Convert the Interial XYZ coordinates to the satellite coordinate system (RSW), then plot the radial component.

Starting Systems: - PCE data - J2000 Coordinate System - Inertial satellite State Vector: \([x, y, z, \dot{x}, \dot{y}, \dot{z}]\) (m) - ORBFIL data - Mean of year Coordinate System - Inertial satellite State Vector: \([x, y, z, \dot{x}, \dot{y}, \dot{z}]\) (m)

Convert from ``XYZ`` to ``RSW``

From Vallado pg. 164:

# %load_ext autoreload
# %autoreload 2

# from PYGEODYNAnalysis_icesat2PCEtrajectory import ARCOVERLAP_2arcs_ObsResids_RSW_radial

# fig = make_subplots(rows=2, cols=1,
#             subplot_titles=(['Radial Component', 'Residual (PCE-ORBFIL)']),
#             vertical_spacing = 0.2,
#             specs=[ [{"secondary_y": False }],
#                     [{"secondary_y": False }]],)

# arc1 = '2018.299'
# arc2 = '2018.299'

# fig = ARCOVERLAP_2arcs_ObsResids_RSW_radial(fig, Obj_Geodyn, 0, arc1, arc2)
# # fig = ARCOVERLAP_2arcs_ObsResids_RSW_radial(fig, Obj_Geodyn3, 1, arc1, arc2)
# # fig = ARCOVERLAP_2arcs_ObsResids_RSW_radial(fig, Obj_Geodyn1, 2, arc1, arc2)

# fig.show(config=config)

# %load_ext autoreload
# %autoreload 2

# from PYGEODYNAnalysis_icesat2PCEtrajectory import ARCOVERLAP_2arcs_ObsResids_NTW_intrack

# fig = make_subplots(rows=2, cols=1,
#             subplot_titles=(['In-Track Component', 'Residual (PCE-ORBFIL)']),
#             vertical_spacing = 0.2,
#             specs=[ [{"secondary_y": False }],
#                     [{"secondary_y": False }]],)

# arc1 = '2018.299'
# arc2 = '2018.299'

# fig = ARCOVERLAP_2arcs_ObsResids_NTW_intrack(fig, Obj_Geodyn, 0, arc1, arc2)
# # fig = ARCOVERLAP_2arcs_ObsResids_NTW_intrack(fig, Obj_Geodyn3, 1, arc1, arc2)
# # fig = ARCOVERLAP_2arcs_ObsResids_NTW_intrack(fig, Obj_Geodyn1, 2, arc1, arc2)

# fig.show(config=config)